Stichprobenverteilung und Schätzung mittels (Re-)Sampling

Prof. Dr. Jörg Schoder

2023-05-29

Urnenmodell

Quelle: moderndive.com




Zufallsexperiment und Stichprobenziehung

Stichprobenziehung

Ergebnis einer Stichprobe

Ergebnis mehrerer Stichproben

Datensatz zum physischen Experiment

library(moderndive)
tactile_prop_red
## # A tibble: 33 × 4
##    group            replicate red_balls prop_red
##    <chr>                <int>     <int>    <dbl>
##  1 Ilyas, Yohan             1        21     0.42
##  2 Morgan, Terrance         2        17     0.34
##  3 Martin, Thomas           3        21     0.42
##  4 Clark, Frank             4        21     0.42
##  5 Riddhi, Karina           5        18     0.36
##  6 Andrew, Tyler            6        19     0.38
##  7 Julia                    7        19     0.38
##  8 Rachel, Lauren           8        11     0.22
##  9 Daniel, Caroline         9        15     0.3 
## 10 Josh, Maeve             10        17     0.34
## # ℹ 23 more rows

Stichprobenverteilung des physischen Experiments

library(tidyverse)
ggplot(tactile_prop_red, aes(x = prop_red)) +
  geom_histogram(binwidth = 0.05, boundary = 0.4,
                 color = "white") +
  scale_y_continuous(limits = c(0, 10), breaks = c(0:10)) +
  labs(x = "Anteil roter Kugeln aus insgesamt 50 (roten und weißen) Kugeln",
       y="Anzahl",
       title = "Verteilung von 33 Anteilswerten roter Kugeln")

Wahre Verteilung

Daten

bowl
## # A tibble: 2,400 × 2
##    ball_ID color
##      <int> <chr>
##  1       1 white
##  2       2 white
##  3       3 white
##  4       4 red  
##  5       5 white
##  6       6 white
##  7       7 red  
##  8       8 white
##  9       9 red  
## 10      10 white
## # ℹ 2,390 more rows

Anzahl und Anteil roter Kugeln

red_true <- bowl %>%
               summarize(Anzahl_rot = sum(color == "red"),
                     Anteil_rot = sum(color == "red")/length(color)
            )
red_true  %>%
    mutate(Anteil_rot=paste0(Anteil_rot*100,"%"))
## # A tibble: 1 × 2
##   Anzahl_rot Anteil_rot
##        <int> <chr>     
## 1        900 37.5%

(Virtuelles) Re-Sampling und Punktschätzung

“Kleine Schaufel” (Stichprobengröße n = 25)

set.seed(23)
n <- 25
rep <- 1000
virtual_samples_25 <- bowl %>%  #Ziehung aus der Urne (Grundgesamtheit)
                        rep_sample_n(size = n,
                                     reps = rep)
virtual_samples_25
## # A tibble: 25,000 × 3
## # Groups:   replicate [1,000]
##    replicate ball_ID color
##        <int>   <int> <chr>
##  1         1     925 red  
##  2         1    2332 white
##  3         1    1272 white
##  4         1    1523 red  
##  5         1     456 white
##  6         1    1707 white
##  7         1    1837 white
##  8         1     913 white
##  9         1    1015 white
## 10         1     149 red  
## # ℹ 24,990 more rows
virtual_prop_red_25 <- virtual_samples_25 %>%
                              group_by(replicate) %>%
                              summarize(red = sum(color == "red")) %>%
                              mutate(prop_red = red / n)
virtual_prop_red_25
## # A tibble: 1,000 × 3
##    replicate   red prop_red
##        <int> <int>    <dbl>
##  1         1    10     0.4 
##  2         2     9     0.36
##  3         3     6     0.24
##  4         4     3     0.12
##  5         5    11     0.44
##  6         6     8     0.32
##  7         7     8     0.32
##  8         8    11     0.44
##  9         9     8     0.32
## 10        10    11     0.44
## # ℹ 990 more rows
virtual_prop_red_25 %>% 
  ggplot(aes(x = prop_red)) +
      geom_histogram(binwidth = 0.05, boundary = 0.4, color = "white") +
      labs(x = paste0("Anteil roter Kugeln (aus ",n,")"),
           title = paste0("Kleine Schaufel (n=",n,")")) +
      geom_vline(xintercept = red_true$Anteil_rot,color='red')

“Mittlere Schaufel” (Stichprobengröße n = 50)

set.seed(23)
n <- 50
virtual_samples_50 <- bowl %>%  # Ziehung aus der Urne (Grundgesamtheit)
                         rep_sample_n(size = n,
                                      reps = rep)
virtual_prop_red_50 <- virtual_samples_50 %>%
                              group_by(replicate) %>%
                              summarize(red = sum(color == "red")) %>%
                              mutate(prop_red = red / n)
virtual_prop_red_50
## # A tibble: 1,000 × 3
##    replicate   red prop_red
##        <int> <int>    <dbl>
##  1         1    20     0.4 
##  2         2    10     0.2 
##  3         3    18     0.36
##  4         4    19     0.38
##  5         5    18     0.36
##  6         6    18     0.36
##  7         7    18     0.36
##  8         8    16     0.32
##  9         9    25     0.5 
## 10        10    19     0.38
## # ℹ 990 more rows
virtual_prop_red_50 %>% 
  ggplot(aes(x = prop_red)) +
      geom_histogram(binwidth = 0.05, boundary = 0.4, color = "white") +
      labs(x = paste0("Anteil roter Kugeln (aus ",n,")"),
           title = paste0("Kleine Schaufel (n=",n,")")) +
      geom_vline(xintercept = red_true$Anteil_rot,color='red')

“Große Schaufel” (Stichprobengröße n = 100)

set.seed(23)
n<-100
virtual_samples_100 <- bowl %>%  # Ziehung aus der Urne (Grundgesamtheit)
                          rep_sample_n(size = n,
                                       reps = rep)
virtual_samples_100
## # A tibble: 100,000 × 3
## # Groups:   replicate [1,000]
##    replicate ball_ID color
##        <int>   <int> <chr>
##  1         1     925 red  
##  2         1    2332 white
##  3         1    1272 white
##  4         1    1523 red  
##  5         1     456 white
##  6         1    1707 white
##  7         1    1837 white
##  8         1     913 white
##  9         1    1015 white
## 10         1     149 red  
## # ℹ 99,990 more rows
virtual_prop_red_100 <- virtual_samples_100 %>%
                              group_by(replicate) %>%
                              summarize(red = sum(color == "red")) %>%
                              mutate(prop_red = red / n)
virtual_prop_red_100
## # A tibble: 1,000 × 3
##    replicate   red prop_red
##        <int> <int>    <dbl>
##  1         1    31     0.31
##  2         2    37     0.37
##  3         3    37     0.37
##  4         4    35     0.35
##  5         5    43     0.43
##  6         6    36     0.36
##  7         7    34     0.34
##  8         8    28     0.28
##  9         9    35     0.35
## 10        10    44     0.44
## # ℹ 990 more rows
virtual_prop_red_100 %>% 
  ggplot(aes(x = prop_red)) +
      geom_histogram(binwidth = 0.05, boundary = 0.4, color = "white") +
      labs(x = paste0("Anteil roter Kugeln (aus ",n,")"),
           title = paste0("Kleine Schaufel (n=",n,")")) +
      geom_vline(xintercept = red_true$Anteil_rot,color='red')

Grundproblem der induktiven Statistik

Bedeutung der Zufallsaufwahl bei der Datenerhebung (!)

Wenn der wahre Wert unbekannt ist

(Virtuelles) Ziehen einer Stichprobe

set.seed(23)
n <- 50
virtual_shovel <- bowl %>% #Ziehen aus der Grundgesamtheit
                     rep_sample_n(size = n)
virtual_shovel
## # A tibble: 50 × 3
## # Groups:   replicate [1]
##    replicate ball_ID color
##        <int>   <int> <chr>
##  1         1     925 red  
##  2         1    2332 white
##  3         1    1272 white
##  4         1    1523 red  
##  5         1     456 white
##  6         1    1707 white
##  7         1    1837 white
##  8         1     913 white
##  9         1    1015 white
## 10         1     149 red  
## # ℹ 40 more rows

Anzahl und Anteil der roten Kugeln

virtual_shovel %>% 
  summarize(num_red = sum(color == "red")) %>% 
  mutate(prop_red = num_red /n)
## # A tibble: 1 × 3
##   replicate num_red prop_red
##       <int>   <int>    <dbl>
## 1         1      20      0.4

Punktschätzung und Erwartungstreue

\[\sigma^2=\frac{p\cdot (1-p)}{n}\]

die vorhandenen Daten wurden als Daten der gesamten Population behandelt, mithin ist der zur Berechnung der Varianz erforderliche wahre Anteilswert (\(p\)) bekannt.

\[\hat{\sigma}^2=\frac{n}{n-1}\sigma^2=\frac{\hat{p}\cdot(1-\hat{p})}{n-1}\]

Statistische Inferenz

Stichprobenverteilung

Reliabilität und Validität

Stichprobenverteilungen…

…und wahrer Wert (rote Linien)

Vergleich der Standardfehler

virtual_prop_red_25 %>%
  summarize(sd = sd(prop_red))
## # A tibble: 1 × 1
##       sd
##    <dbl>
## 1 0.0935
virtual_prop_red_50 %>%
  summarize(sd = sd(prop_red))
## # A tibble: 1 × 1
##       sd
##    <dbl>
## 1 0.0673
virtual_prop_red_100 %>%
  summarize(sd = sd(prop_red))
## # A tibble: 1 × 1
##       sd
##    <dbl>
## 1 0.0486

Zentraler Grenzwertsatz

Mit zunehmendem Stichprobenumfang nähert sich (1) die Stichprobenverteilung eines Punktschätzers (bspw. Anteil roter Kugeln) einer Normalverteilung an und nimmt (2) die Streuung dieser Stichprobenverteilungen ab ( kleinere Standardfehler).

Die Normalverteilung der Stichprobenverteilung resultiert dabei unabhängig von der Verteilung in der Grundgesamtheit - solange die für die Stichproben jeweils berechneten Mittelwerte auf unabhängigen Zufallsvariablen basieren. Die Unabhängigkeit wird dabei durch die Zufälligkeit der Stichprobenziehung sichergestellt.

Beispiele, vgl. Skript Wahrscheinlichkeitstheorie & Zufallsvariablen

Mögliche Punktschätzer

Im bisher verwendeten Beispiel ging es um den Anteil roter Kugeln, mithin um die Schätzung eines Anteilswerts. Das Konzept der Stichprobenziehung kann jedoch auch zur Schätzung anderer Parameter einer Grundgesamtheit verwendet werden:

Mögliche Punktschätzer auf Basis von Stichproben
Fall Parameter der Grundgesamtheit Notation Punktschätzung Symbol(e)
1 Anteil in Grundgesamtheit \(p\) Anteil in Stichprobe \(\widehat{p}\)
2 Mittelwert der Grundgesamtheit \(\mu\) Stichprobenmittelwert \(\overline{x}\) oder \(\widehat{\mu}\)
3 Differenz von Anteilen einer Grundgesamtheit \(p_1 - p_2\) Differenz von Anteilen einer Stichprobe \(\widehat{p}_1 - \widehat{p}_2\)
4 Differenz von Mittelwerten der Grundgesamtheit \(\mu_1 - \mu_2\) Differenz von Stichprobenmittelwerten \(\overline{x}_1 - \overline{x}_2\) oder \(\widehat{\mu}_1 - \widehat{\mu}_2\)
5 Empirischer Regressionskoeffizient (Grundgesamtheit) \(\beta_1\) Angepasster (“fitted”) Regressionskoeffizient (Stichprobe) \(b_1\) oder \(\widehat{\beta}_1\)




Bootstrapping und Konfidenzintervalle

Bootstrapping und Intervallschätzungen

Stichproben-Variation ohne Kenntnis der Grundgesamtheit

Bootstrapping als Methode zur Untersuchung der Stichproben-Variation auf Basis einer (einzigen) Stichprobe.

Punktschätzung vs. Intervallschätzung

Bootstrapping als Methode zur Ermittlung von Konfidenzintervallen ohne die Notwendigkeit zur Annahme einer bestimmten Verteilung der Grundgesamtheit.

Stichprobenverteilung (Reminder)

set.seed(23)
virtual_samples <- bowl %>% # Ziehen aus der Grundgesamtheit
                      rep_sample_n(size = 50, reps = 1000)
sampling_distribution <- virtual_samples %>%
                                  group_by(replicate) %>%
                                  summarize(red = sum(color == "red")) %>%
                                  mutate(prop_red = red / 50)

## # A tibble: 1 × 1
##       se
##    <dbl>
## 1 0.0673

Bootstrapping-Verteilung

Einzelne Stichprobe aus dem physischen Experiment

bowl_sample_1
## # A tibble: 50 × 1
##    color
##    <chr>
##  1 white
##  2 white
##  3 red  
##  4 red  
##  5 white
##  6 white
##  7 red  
##  8 white
##  9 white
## 10 white
## # ℹ 40 more rows
stats_sample_1 <- bowl_sample_1 %>%
                  summarize(Anzahl_rot=sum(color=='red'),
                            Anteil_rot=sum(color=='red')/
                                                length(color))
stats_sample_1
## # A tibble: 1 × 2
##   Anzahl_rot Anteil_rot
##        <int>      <dbl>
## 1         21       0.42

In der Stichprobe von Ilyas und Yohan befinden sich insgesamt 21 rote Kugeln, d.h. der Anteil roter Kugeln entspricht in ihrer Stichprobe 42%.

Bootstrapping-Verteilung

set.seed(23)
bowl_sample_1 %>%  #Stichprobendaten von Ilyas und Yohan
  rep_sample_n(size = 50, reps = 1000) %>% 
    group_by(replicate) %>%
    summarize(red = sum(color == "red")) %>%
    mutate(prop_red = red / 50)
## # A tibble: 1,000 × 3
##    replicate   red prop_red
##        <int> <int>    <dbl>
##  1         1    21     0.42
##  2         2    21     0.42
##  3         3    21     0.42
##  4         4    21     0.42
##  5         5    21     0.42
##  6         6    21     0.42
##  7         7    21     0.42
##  8         8    21     0.42
##  9         9    21     0.42
## 10        10    21     0.42
## # ℹ 990 more rows
set.seed(23)
bowl_sample_1 %>%  #Stichprobendaten von Ilyas und Yohan
  rep_sample_n(size = 50, reps = 1000, replace=TRUE) %>%  # Ziehen mit Zurücklegen
    group_by(replicate) %>%
    summarize(red = sum(color == "red")) %>%
    mutate(prop_red = red / 50)
## # A tibble: 1,000 × 3
##    replicate   red prop_red
##        <int> <int>    <dbl>
##  1         1    20     0.4 
##  2         2    22     0.44
##  3         3    21     0.42
##  4         4    21     0.42
##  5         5    21     0.42
##  6         6    24     0.48
##  7         7    24     0.48
##  8         8    24     0.48
##  9         9    21     0.42
## 10        10    25     0.5 
## # ℹ 990 more rows

Bootstrapping mit dem infer-Paket

library(infer)
set.seed(23)
bootstrap_distribution <- bowl_sample_1 %>%
                              specify(response = color,
                                      success = "red") %>%
                              generate(reps = 1000,
                                       type = "bootstrap") %>%
                              calculate(stat = "prop")
bootstrap_distribution %>%
                  visualize() +
                  labs(x = "Anteil roter Kugeln (aus 50)",
                       y= "Anzahl",
                       title = "Bootstrapping-Verteilung") +
                  geom_vline(xintercept = stats_sample_1$Anteil_rot,
                             color='red', linetype='dashed')

## # A tibble: 1 × 1
##       se
##    <dbl>
## 1 0.0701

Ermittlung von Konfidenzintervallen

Nutzung der Funktionen im infer-Paket

Schritt 1: specify()

#bowl_sample_1 %>%        
#    specify(response = color)   # funktioniert nicht - "success" (also das "Ereignis A") muss definiert werden!

bowl_sample_1 %>%
    specify(response = color, success = "red")
## Response: color (factor)
## # A tibble: 50 × 1
##    color
##    <fct>
##  1 white
##  2 white
##  3 red  
##  4 red  
##  5 white
##  6 white
##  7 red  
##  8 white
##  9 white
## 10 white
## # ℹ 40 more rows

Schritt 2: generate()

bowl_sample_1 %>%
  specify(response = color, success = "red") %>%
  generate(reps = 1000, type = "bootstrap")
## Response: color (factor)
## # A tibble: 50,000 × 2
## # Groups:   replicate [1,000]
##    replicate color
##        <int> <fct>
##  1         1 white
##  2         1 white
##  3         1 white
##  4         1 red  
##  5         1 white
##  6         1 white
##  7         1 red  
##  8         1 red  
##  9         1 white
## 10         1 white
## # ℹ 49,990 more rows

Schritt 3: calculate()

sample_1_bootstrap <- bowl_sample_1 %>%
                           specify(response = color,
                                   success = "red") %>%
                           generate(reps = 1000,
                                    type = "bootstrap") %>%
                           calculate(stat = "prop")
sample_1_bootstrap
## Response: color (factor)
## # A tibble: 1,000 × 2
##    replicate  stat
##        <int> <dbl>
##  1         1  0.48
##  2         2  0.48
##  3         3  0.4 
##  4         4  0.4 
##  5         5  0.42
##  6         6  0.54
##  7         7  0.5 
##  8         8  0.4 
##  9         9  0.42
## 10        10  0.48
## # ℹ 990 more rows

Schritt 4: visualize()

## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1      0.3     0.54
sample_1_bootstrap %>%
        visualize(bins = 15) +
        shade_confidence_interval(endpoints = percentile_ci_1) +
        geom_vline(xintercept = 0.42, linetype = "dashed")

Interpretation Konfidenzintervall

“Mit 95%iger Wahrscheinlichkeit enthält das Konfidenzintervall den wahren Wert \(p\).”

“Bei hinreichend großer Anzahl wiederholter Stichprobenziehungen, können wir erwarten, dass etwa 95 % der resultierenden Konfidenzintervalle den wahren Wert des Populationsparameters umfassen/beinhalten.”

Perzentil-Methode

Standardfehler-Methode

\[\hat{p} \pm 1.96 \cdot SE = (\hat{p} - 1.96 \cdot SE, \, \hat{p} + 1.96 \cdot SE)\] * Der Standardfehler kann dabei… * …mittels Bootstrapping ermittel werrden (s.o.) * …oder basierend auf der theoretischen Verteilung:

\[\text{SE}_{\widehat{p}} \approx \sqrt{\frac{\widehat{p}(1-\widehat{p})}{n}}\]